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Abstract 

Background: Translation efficiency of certain mRNAs can be regulated through a cytoplasmic polyadenylation 
process at the pre-initiation phase. A translational regulator controls the polyadenylation process and this 
regulation depends on its posttranslational modifications e.g., phosphorylation. The cytoplasmic polyadenylation 
binding protein (CPEB1) is one such translational regulator, which regulates the translation of some mRNAs by 
binding to the cytoplasmic polyadenylation element (CPE). The cytoplasmic polyadenylation process can be turned 
on or off by the phosphorylation or dephosphorylation state of CPEB1. A specific example could be the regulation 
of Calcium/Calmodulin-dependent protein kinase II (aCaMKII) translation through the phosphorylation/ 
dephosphorylation cycle of CPEB1. 

Result: Here, we show that CPEB1 mediated polyadenylation of aCaMKII mRNA can result in a bistable switching 
mechanism. The switch for regulating the polyadenylation is based on a two state model of aCaMKII and its 
interaction with CPEB1. Based on elementary biochemical kinetics a high dimensional system of non-linear ordinary 
differential equations can describe the dynamic characteristics of the polyadenylation loop. Here, we simplified this 
high-dimensional system into approximate lower dimension system that can provide the understanding of 
dynamics and fixed points of original system. These simplified equations can be used to develop analytical 
bifurcation diagrams without the use of complex numerical tracking algorithm, and can further give us intuition 
about the parameter dependence of bistability in this system. 

Conclusion: This study provides a systematic method to simplify, approximate and analyze a translation/activation 
based positive feedback loop. This work shows how to extract low dimensional systems that can be used to obtain 
analytical solutions for the fixed points of the system and to describe the dynamics of the system. The methods 
used here have general applicability to the formulation and analysis of many molecular networks. 



Background 

Cellular signaling pathways that can operate in a switch 
like manner are called bistable systems [1,2]. A bistable 
system has the ability to switch between two distinct 
stable steady states and such a system cannot rest in 
any intermediate state [3,4]. In response, to an external 
stimulus a bistable system can move from one state to 
another. If this switching is permanent then such a system 
is called irreversible otherwise it is a reversible switch [5]. 
The bistability in a signaling network is typically due to a 
positive feedback loop or double negative feedback loop 
[3]. However, the presence of a positive or a double nega- 
tive feedback loop does not guarantee bistability [4]. 
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In addition, to these feedback loops a biological network 
must have non-linear interactions to exhibit a bistable 
behavior. Previous, experimental work has described sev- 
eral examples of naturally occurring bistable system 
[6-18]. Still bistability is not considered to be a unifying 
theme of cellular signaling networks and more experimen- 
tal work is needed to establish bistability as one of the gen- 
eral mechanism of cell signaling [4]. Typically, bistable 
biological systems were described either at the level of 
gene expression due to the regulation of gene expression 
by transcription factors, or at the level of posttranslational 
modifications e.g., activation-deactivation cycle due to 
phosphorylation. Here, we present a model of bistability 
that can arise from the control of gene expression at the 
level of translation of new proteins 
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Cytoplasmic polyadenylation regulates the translation 
efficiency of certain mRNAs through modulating the 
length of 3' poly (A) tail [19,20]. Polyadenylation can be 
regulated through a short nucleotide sequence known as 
cytoplasmic polyadenylation element (CPE) located in 
their 3' UTR. A CPE binding protein (CPEB1), represses 
the translation through its dual interactions with CPEs 
and other mRNA binding proteins. Phosphorylation of 
CPEB1 changes its interactions with these other pro- 
teins, promotes polyadenylation and increases the length 
of the poly (A) tail [19,20]. The mRNAs with a longer 
poly (A) tail are more likely to be translated compared 
to mRNAs with a shorter poly (A) tail [21,22]. The 
exact mechanism through which the longer poly (A) tail 
enhances the translation efficiency is not clear. However, 
it is believed that longer poly(A) tail leads to a circular 
mRNA which enhances the translation efficiency 
through recycling the translation machinery on the 
mRNA frame thus increasing the possibility of mRNA 
translation through initiation [23,24]. 

The translation of a highly abundant brain protein 
aCaMKII (2% of total brain protein is aCaMKII) is regu- 
lated through activity induced polyadenylation. The 
aCaMKII-mRNA contains two CPE elements in its 
3'UTR and its translation can be regulated through phos- 
phorylation of CPEB1 [25,26]. Recent studies have shown 
that aCaMKII can phosphorylate CPEB1 and therefore, 
possibly modulate its own translation through a positive 
feedback loop. Here, we examine the hypothesis that the 
positive feedback loop between aCaMKII and CPEB1 
forms a bistable switch which regulates the translation of 
aCaMKII. The aim of this paper is to obtain analytical 
expressions for the bistability of the CPEB1- aCaMKII 
molecular pair as a generic example for such systems. 

In this paper, we analyze the mathematical properties of 
this molecular loop. The characteristics of this loop are 
analyzed by evaluating the dynamics and directly locating 
the fixed points. Using the elementary biochemical 
kinetics we develop a molecular model of self-sustained 
polyadenylation based translation loop. This simple mole- 
cular model is represented by six differential equations. 
We simplify this model through introducing an approxi- 
mation and algebraic manipulations. These, systematic 
simplifications result in a three dimensional differential 
equation based model. Further approximations can reduce 
this to a single dynamical equation. Based on our approxi- 
mate equation we also developed an approximate analyti- 
cal method to directly locate the fixed points of this 
system. We compare these analytical results to the numer- 
ical bifurcation diagrams obtained through numerical 
tracking of the complete system of equations and show 
their correspondence. Our results demonstrate that such a 
positive feedback loop which involves the control of 



translation through polyadenylation can indeed be bistable 
over a wide range of parameters. This simplified model, 
though motivated by the aCaMKII-CPEBl loop, could be 
seen as a generic model for such a positive feedback loops 
which involves translation, and degradation of proteins. 
Such feedback loops do not conserve the quantity of these 
proteins and are therefore qualitatively different than most 
post translational feedback models [3,4]. Since, we use a 
simplified model we can obtain approximate analytical 
results, which provide us with intuition about how such 
feedback systems operate. 

Method 

A. Complete Model Equations 

The following set of reactions is used to describe the 
interactions between the aCaMKII and CPEB1 mole- 
cule. These biochemical reactions are based on elemen- 
tary Michalis-Menten type kinetics. The dynamical 
variable X represents aCaMKII and Y represents the 
CPEB1. The P subscript represents the phosphorylated 
form and an A as superscript represents the active form. 



x+(Ca +2 ) 4 .CaM 



x P 



X + Xp 

Xp + P 
Y + xp 

Yp + P 
Y P +T 



Ci 



2xp 



X + P 



Yp+Xp 



k 7 



Y + P 



k8 



c. 



Yp + X + T 



(Rl) 



(R2) 



(R3) 



(R4) 



(R5) 



(R6) 



From above reactions following differential equations 
can be deduced. 



dX P 



^ - -ki * X * X P + k 2 * Ci + 2 * k 3 * Ci — k 4 * Xp * P + k 55 * C2 + k 6 * C2 
-k 5 * Y * X P + k 10 * x * U - k 10 i 0 * X P — A.! * (X P - Xp ba sal) 



dY P 

~dT : 



-k 7 * Yp * P + k6 * C2 — ks * Yp * T + kss * C3 + kg * C3 (A — 2) 



-X 2 * (X - X basa i) 



°* Xp (A -3) 
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dCi 

— — = ki * X * Xp — k 2 * Ci — k3 * Ci 
dt 



dC 2 

— — = k 3 * Y * Xp — kss * C 2 — kg * C2 
dt 



dC 3 
~dT 



kg * Yp * T — k 88 * C 3 — k 9 * C 3 



(A -4) 
(A -5) 
(A -6) 



B. Model analysis and reduction 

In the following equations, the dynamical variable X 
represents aCaMKII and Y represents the CPEB1. The 
P subscript represents the phosphorylated and active 
form. By using the pseudo-steady state assumptions the 
differential equations representing the complexes [Q- 
C3] can be eliminated: 



Ci 



C 3 = 



ki * X * Xp 
(k 2 + k 3 ) 



ks * Y * Xp 
(k 55 + k 6 ) 



k 8 * Y P * T 
(k 88 + k 9 ) 



(B-l) 



(B-2) 



(B-3) 



The unphosphorylated CPEB1 (Y) is related to phos- 
phorylated CPEB1 (Y P ). 



Y = Y T - Y P 



(B-4) 



Where, Y T is the total amount of CPEB1. The aCaM- 
KII molecules are either in free or in bound form, there- 
fore, the total concentration of aCaMKII is given by 
following equation. 



X T = X + X P + 2 * Ci + C 2 



(B-5) 



The differential equation representing the phosphory- 
lated CPEB1 (Y P ) from Rj-Rg 



dY P 



-k 7 * Yp * P + k 6 * C2 — k 8 * Yp * T + k 88 * C3 + k 9 * C3 (B — 6) 

By requiring a steady state ^Xli = 0 and substituting 

dt 

k6*k5 
a = — — r we get 



(k 55 + k 6 ) 

a * Yt * Xp 



Y P 



(k 7 * P + a * Xp) 



(B-7) 



Also from (B-l, B-2, B-3, B-4, B-5 and B-7) we can 
obtain the value of X. 



X T - X P - N * X P * Y T 



X = 



a * Xp * Yp \ 
k 7 * P + a * X P ) (B - 8) 



(1 + 2*M*X P ) 



Where, N and M are constants defined as 

k 5 ki 

N = ,M = 

(k 55 +k 6 ) (k 2 +k 3 ) 

As described in the result section the equation 2 is as 

follows: 



H(C 3 ,X T ) = k 9 *C 3 -A*X T 
Substituting C3 from (B-3) in [2] 
k 8 * T 



H(Y P ,X T ) = k 9 * 



*Y P 



(2) 



■A.*X T (B-9) 



k 88 + k 9/ 
Placing Y P from (B-7) into (B-9) 

H (x P/ x T ) = k 9 *f^)*f^|^)-x*x T (B-10) 
Vk 88 + k 9 y V a * x p + p 7/ 



Where, P 7 = P*k 7 and P 4 = P*k 4 . Defining, U = (Ca 
k 9 * k 8 * T 



+2\ 



4-CaM, 
c = Y T * 



(k 88 + k 9 ) 
k 9 * k 8 * T 



and 



= b*Y 7 



or 



(k 88 + k 9 ) 
We obtain 

h(Xp ' Xt) = (x^)- a * Xt (B " n) 

The equation 1 describes the rate of change of total 
concentration of aCaMKII (X T ) is obtained from the 
balance between the new synthesis of aCaMKII and its 
degradation. This equation is further transformed to the 
approximate solution as shown by equation 4. Similarly, 
from R1-R6 an equation representing the dynamics of 
phosphrylated aCaMKII can be constructed. 



dX P 



= — ki * X * Xp + 1<2 * Ci + 2 * k3 * Ci — 1<4 * Xp * P + kss * C2 + kg * C2 
-k 5 *Y*Xp + k 1 o*X*U-k 1 oio*X P 



(B - 12) 



The equation B-12 is further simplified in only two 



variables i.e., X T and X P by placing the 



dY P 
~dT 



0 and put- 



ting the values of X, Q, C 2 , Y, and Y P from B-l to B-8 
and further simplifying we get the following expression. 

X T = Ii + h + h + U (B - 13) 

The I1J2J3 and I 4 are defined as follows: 

Ii = Di. *Xp./(ai.*X P + a 2 ) (B - 14) 

h = [Di.*X P . 2 ./(k 3 .*X P ./2 + k 3 .*a 2 ./(2*ai))](B- 15) 



h = Xp 



(B - 16) 
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I 4 =Xp.*Y T ./(k 6 ./d 1+ Xp./k 7 ) 



(B - 17) 



The expression in B-14 to B-17 are defined as k 7 = k 7 
•P, k 4 = k 4 *P, ai = (k 3 % )./(k 2 + k 3 ) a 2 = kio *U, D l = 
k 4 + kioio , di = (k 6 .*k 5 )./(k 55 + k 6 ) 

The expression of B-13 gives us a function where X T 
= f(X P ), It is almost impossible to invert this function to 
get X P = f _1 (X T )so we numerically approximated it 
through function fitting. 

Equation B-12 with the aid of equation Bl-B 13 can 
be transformed to following polynomial in terms of X P . 



(Xp) 4 + Z12 * (Xp) 3 + Z13 * (Xp) 2 + Z14 * Xp = 0 


(B- 


18) 


zi = a * c * M * k 3 * P 7 


(B- 


19) 


z 2 = a * c * kio * U * P 7 


(B- 


20) 


z 3 = a 2 * c * M * k 3 


(B- 


21) 


z 4 = a 2 * c * F * kio * U 


(B- 


22) 


Z5 = 2*P4*M*F + 2*M* kioio *F + F*M*k3 


(B- 


23) 


z 6 = 2*P4*P 7 *M+2*M*P7*kioio+F*P4+F*k 10 io+M*k 3 *P7+F*kio*U+N*M*k3*YT 


(B- 


24) 


z 7 = P 7 * P 4 + P 7 * kioio + kio * U * P7 + N * kio * U * P7 * Yt (B — 


25) 


z§ = a * Z5 


(B- 


26) 


z 9 = a * z 6 + P 7 * Z5 — z 3 


(B- 


27) 


Z10 = a * z 7 + P 7 * Z6 — zi — z 4 


(B- 


28) 


zn = p 7 * z 7 - z 2 


(B- 


29) 



Z12 



Zl3 



Z9 
Z 8 

zip 
z 8 



zn 

Zl4 = 

z 8 



(B - 30) 



(B - 31) 



(B - 32) 



Where, the coefficients z 12 >Zi3 and z 14 are defined by 
following expressions (B-19-B-32). The equation (B-18) 
is obtained through further simplifying (B-13) such that 
a fourth order polynomial is obtained in terms of X P . 
This 4'th order polynomial has an analytical solution 
because it has no zero order term and therefore has one 



solution X P = 0, and the other solutions are the solu- 
tions of a third order polynomial 

Results 

Analysis of the model of aCaMKII synthesis through 
polyadenylation 

Our simplified model of a self-sustained polyadenylation 
of aCaMKII-mRNA (Figure 1) is based on biochemical 
interactions between a plasticity related kinase aCaMKII 
and its translational regulator CPEB1 through a positive 
feedback. This model of polyadenylation based transla- 
tion of aCaMKII-mRNA (Figure 1) is composed of two 
molecular components which interact through a closed 
loop. 1) The aCaMKII protein which can be in two 
states: inactive, and phosphorylated, active. 2) The 
CPEB1 a translational regulator to regulate the polyade- 
nylation, can be either in the phosphorylated or unpho- 
sphorylated state. The phosphorylated CPEB1 promotes 
the translation at pre-initiation phase through polyade- 
nylation. Here, in this simple model we assume that 
CPEB1 is phosphorylated only by active and phosphory- 
lated aCaMKII. In this molecular scheme the aCaMKII 
protein is removed at a certain degradation rate [27]. 
The synthesis of new aCaMKII protein regulated by 
polyadenylation provides the necessary compensation 
for the amount removed due to degradation. Our model 
shows how the concentration of aCaMKII can be main- 
tained at multiple levels despite the synthesis of new 
molecules and removal due to protein degradation. 

This system is described by a set of differential equa- 
tion (Method A). We use notation in which CaMKII is 
denoted as X, and phosphorylated aCaMKII as X P . 

Similarly CPEB is denoted as Y and phosphorylated 
CPEB1 as Y P . On purpose we have chosen a simplified 
version of these components, when compared to our 
previous work [28] in order to gain intuition into the 
behavior of the system while keeping the bistability. We 
can easily show that for appropriate parameters this sys- 
tem is bistable. However, gaining an intuitive under- 
standing of this bistability is difficult even in this 
simplified system, because the dimensionality is still too 
large. 

The two-state assumption of the aCaMKII is a signifi- 
cant simplification compared to more realistic multi- 
state models and was chosen to enable simpler analysis. 
The synthesis of new proteins in our proposed molecular 
model is controlled by the phosphorylated CPEB1 mole- 
cules. The translation model is also a simplified represen- 
tation of a complex molecular system. The aim of this 
paper is to generate a much reduced dynamical model 
that can nevertheless capture the qualitative behavior of 
the system, and provide us with a more intuitive under- 
standing of dynamical behavior of the system. 
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(Ca +2 ) 4 ,CaM c 




Figure 1 The simplified model of CPEB1 mediated polyadenylation of aCaMKII through a self-sustaining ocCaMKII-CPEB 1 molecular 
loop. The aCaMKII, molecule can be inactive, or active and phosphorylated state, while CPEB1 can be in unphosphorylated and phosphorylated 
states. The aCaMKII molecule in both these states degrades with a certain degradation rate. The active and phosphorylated CaMKII 
phosphorylate the CPEB1 molecule, which through polyadenylation of aCaMKII mRNA generates the new CaMKII molecule. Here, the synthesis 
of new aCaMKII molecule is represented by a translation step (T) which, describe the multi-step translation process through a single biochemical 
reaction R6 as described in appendix. 



One level of simplification is to set the derivatives of 
the three complexes (Q-C3) to zero (Method B), to 
obtain a pseudo steady state approximation, similar to 
that used in the standard michaelis-menten approxima- 
tion. This reduces our system to 3 dynamical equations, 
and three functional expressions, as described by equa- 
tions I, B-6, and B-12. This approximation simplifies 
the dynamics, and can be used to more easily find the 
steady states of this system numerically. We can further 
approximate this system by a single differential equation. 
The key simplifying assumption required here is the fol- 
lowing equation: 

= k 9 * C3 — A * Xt (1) 

dt 

Where, C 3 , as described by equation B-3 in the 
method B, is the concentration of the phosphorylated 
form of CPEB1, bound to the translation machinery, or 
in other words the concentration of active translation 
machinery available for producing new aCaMKII. Here 
k 9 is the forward rate of generating new aCaMKII. This 
equation assumes that aCaMKII degrades at the same 
rates both in free and bound forms. The X T in above 
equation (Eq. 1) is the concentration of total aCaMKII. 
In our detailed system of equations (Method A) we do 
not assume the degradation of complexes. However, this 



approximation holds if the relative concentration of 
complexes at steady state is small. 
dX T 

Defining = H(C3,Xx) we have that: 

dt 

H(C 3 , X T ) = k 9 * C 3 - k * X T (2a) 

By using a pseudo steady state approximation on some 
of the faster dynamical variables, we obtain (Method B) 
that: 

H(X P/ X T ) = ( C * Xp - 1 - A*X T (3) 

Where the parameters ( c\ 'P7' and V, are defined in 
Method B. The parameter c is linearly dependent on Y T , 
and is inversely proportional to degradation rate A. It 
shows that amount of total aCaMKII increases with 
increase in translation (either due to increase in transla- 
tion rate or concentration of translation machinery "T") 
and decreases with increase in degradation rate. The 
Kmi/2 of this process is P 7 /a and depends on phos- 
phatses and dephosphorylation rate of CPEB1. 

We call the first function of RHS of equation 3 the 
synthesis function and the second one a degradation 
function. We can rewrite equation 3 as 

H(X P/ X T ) + G(X P ) - F(X T ) (4) 
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Where, G is a synthesis function, and F a degradation 
function, and the dynamics are then understood as a 
balance between synthesis and degradation. If we can 
derive a function that relates X p to X T , i.e., X P = g (X T ), 
then equation 4 could be converted to equation of a sin- 
gle dynamical variable: 

H(X T ) = G(g(X T )) - F(X T ) = G'(X T ) - F(X T ) (5) 

and we would have reduced the 6 coupled ODE's to a 
single ODE: 

dX T 

= G'(X T ) - F(X T ) (6) 

at 

Although we have an exact expression for X T as a 
function of X P [(Method B equation 13)], it is difficult 
to analytically invert it to obtain an analytical expression 
for X P = g ( X T ). Instead it is easy to invert this numeri- 
cally using a higher order polynomial fit. 

Using this fit. 

9 

X P = ^a i *X i T = h(X T ) (7) 

i=l 

The fitting performance of this equation is shown in 
additional file 1, Figure Si. The solid black line in addi- 
tional file 1, Figure S x represents the original function 
{Method B B-13}, while the dotted red line represents 
the fit from equation 7. Thus we reduce the 6 differen- 
tial equations to a single approximate ODE of the form 
described in equation 6, where: 

\h(X T ) + P 7 /a/ 

Here, the X P is obtained from fitting function 7 and F 
(X T ) is given by following equation. 

F(X T ) = A * Xx (9) 

In order, to gain some intuitive understanding of this 
system we graphically plot the two terms in equation 6 
(Figure 2a). The function G is shown as dotted red line 
and F as solid blue line. This graphical representation 
(Figure 2a) provides a very simple intuitive explanation of 
the behavior of the system. The function G' can be seen 
as the synthesis function (a source) and the function F as 
the degradation function (a sink). The intersections 
between these two curves are the systems fixed points. 
For low values of X T there is hardly any synthesis, and 
degradation dominates, so the system converges to the 
low fixed point at zero. At higher values of X T there is an 
abrupt rise in the level of protein synthesis, which quickly 
saturates. The cross between G and F in the quickly ris- 
ing portion of G is the unstable fixed point, below it 
degradation dominates and above it synthesis dominates, 



and the level of X T increases until it reaches the third 
intersection, which is the upper stable fixed point. Due to 
the saturation of G, at higher levels of X T , degradation 
dominates again, resulting in a convergence from above 
to the upper stable fixed point. Changing the degradation 
rate will simply change the slope of the F function, and 
changing synthesis and activation parameters, will quan- 
titatively change the shape of the G' function. 

Next we show the impact of changing some of the sys- 
tem parameters on the steady state solution characteristics 
of equation 6 (Figure 2b and 2c). Changing the degrada- 
tion rate only affects the degradation curve (Figure 2b 
solid blue line). As the degradation rate is increased from 
slower (Figure 2b solid blue line 1) to much faster (Figure 
2b solid blue line 4) the system moves from robustly bis- 
table to a mono-stable system. Thus, at much faster degra- 
dation the amount of CaMKII generated is not enough to 
compensate for the loss due to large degradation rate 
therefore, system has only lower stable steady state solu- 
tion (Figure 2b solid blue line 4) as the CaMKII degrada- 
tion rate is decreased a balance between new synthesis 
and degradation is restored and system becomes a bistable 
switch (Figure 2b solid blue lines 2-4). As the degradation 
rate is decreased, and the system as shifts from mono- 
stable to bistable, the lower fixed point does not change, 
the value of X T at the unstable fixed point is increased, 
and the value of X T at the upper fixed point first rapidly 
increases and then plateaus (Figure 2b solid blue line 
2,3,4). 

Other parameters affect only the synthesis curve (G')« 
One such example is parameter k 5 , which quantifies the 
CaMKII mediated activation of CPEB1 (Figure 2c). 
When this activation parameter is set at low value (Fig- 
ure 2c dotted red line 4) the system has only lower 
stable steady state solution. However, as the activation 
rate of CPEB1 is increased, more active CPEB1 is avail- 
able, which in turn stimulates new protein synthesis 
thus shifting the synthesis curve in upward direction 
(Figure 2c dotted red line 3,2,1) and shifting the sys- 
tem from mono-stable to robustly bistable state. 

If we are only interested in the steady states of this 
system we can set all the derivatives to zero and obtain 
the following forth order polynomial to describe the 
steady states of the system: 

(X P ) 4 + zi 2 * (X P ) 3 + zis * (X P ) 2 + zu * X P = 0 (10) 

Where, the coefficients z 12 -z 14 are defined in method B. 
The derivation leading to this result is explained in 
method B. This equation is more precise than equation 6 
because here we avoided the need to use the polynomial 
fit for the inversion. However, this equation can only 
account for the fixed points, and is unable to account for 
system dynamics. This 4'th order polynomial shows that 
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x T 



Figure 2 The approximate analytical solution of ocCaMKII-CPEBI molecular loop. It is obtained through graphically locating all the steady 
state solutions of single differential equation 1. This equation describes the rate of change of total concentration of aCaMKII (X T ) and its 
functional relation between the new synthesis and degradation of CaMKII. The steady state form of equation 1 is transformed into two functions 
F (X T ) and G (X T ). The function F (X T ) (solid blue curve) basically represents the net aCaMKII degradation, whereas the modified function 
G'(X T ) = G [g (X T )] (dotted red curve) represent the net generation of aCaMKII. (a) This analytical solution is developed by setting all parameters 
as described in table 1. Graphs of both functions intersect at three locations which are characterized as steady state solutions of equation 1 and 
provide the approximate analytical solution. The upper steady state solution is located at X T = 95, while the unstable solution is located at X T = 
9.4, and lower steady state solution at X T = 0.0001. (b) The effect of degradation rate on analytical solution. Four different degradation rates are 
selected X = 0.00006 s- (solid blue line #1, bistable system), X, = 0.0001 s 1 (solid blue line #2 bistable system), X = 0 0003 s (solid blue 
line #3 bistable system) and X = 0.0006 s" 1 (solid blue line #4 mono-stable system), (c) The effect of CPEB1 activation parameter k 5 on 
analytical solution. Four different values of CPEB1 activation are selected k 5 = 0.00001 ^/M" 1 .s" 1 (dotted red line #1, mono-stable system), k 5 = 
O.OOOOe^fvrV 1 (dotted red line #2 bistable system), k 5 = 0.0001 \M\s A (dotted red line #3 bistable system) and X = 0.0072 ^M"Y 1 (dotted 
red line #4 bistable system). 



this system has 4 steady states, one is always at zero 
because there is no zero order term. The signs of the 
coefficients determine how many real solutions this poly- 
nomial could have [29] and biologically we are only inter- 
ested in positive real solutions. 

Comparison of the different approximations to the full 
model 

Our aim here is to describe the behavior of the synthesis- 
degradation loop (Figure 1) and the various approxima- 
tions used to simplify it. We can integrate the dynamics 
either using the full system of ODE's (Method A, A1-A6) 
or by using the reduced three dimensional systems 
(Method B, equations 6, B-6, B12), or the approximate 
one dimensional system (equation 6). The fixed points, 
and their parameter dependence can be found either by 
numerical bifurcation analysis of the full system (Method 
A, A1-A6), the reduced three dimensional system (6, B-6, 



B12), the approximate ID system (equation 6) or by the 
4'th order polynomial [10]. The first three methods can be 
used for obtaining the system dynamics and the last 
method only for the fixed points. 

In order, to further compare the three dynamical equa- 
tion methods we first analyzed the dynamical properties 
of this system and then the steady state solution charac- 
teristics of this loop with respect to certain system para- 
meters. The dynamics of this system (Figure 3) is 
analyzed with two different set of parameters. The first 
set of parameters (Condition I) represent a case where 
the approximation of equation 1 is a good aproximation 
(Figure 3a), whereas the second set of parameters (Condi- 
tion II) represents a case where this approximation does 
not hold (Figure 3b). For condition I the approximation 
holds since the amount of aCaMKII bound in biochem- 
ical complexes is negligible, whereas for condition II a 
significant portion is trapped in biochemical complexes. 
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Figure 3 The dynamic characteristics of aCaMKII-CPEB 1 molecular loop. These characteristics are developed through three alternative 
models. First model is based on differential equations A1-A6 (Full model, represented by solid red line (up state) and solid blue line 
(down state)). Second model is based on differential equation 1, B-6 and B-9 (Reduced model, represented by dotted red line (up state) 
and dotted blue line (down state)). Third model is based only on a single differential equation (equation 6) (Single differential equation 
model, represented by dotted black line (up state) and dotted blue line (down state)). System moves from down to up-regulated state 
through a 10 second application of (Ca +2 ) 4 -CaM pulse. The dynamic simulations are carried out with two different set of parameters. For one 
set of parameters the dynamic from three alternative models converge to same upper and lower steady states (a). In this case the dynamics 
converge to same upper steady state at X T = 95 and lower steady state at X T = 0.0 For second selection of parameters the dynamic results from 
three alternative models does not converge to same upper steady state (b) although the lower steady state remains the same for all three 
models. Here, the dynamics from full model converge to upper steady state at X T = 360 and lower state at X T = 0.0001 (Figure 3 b, solid red 
line (up state) and solid blue line (down state)), whereas the dynamics from reduced model (dotted red line (up state) and dotted blue 
line (down state)) and single differential equation model (dotted black line (up state) and dotted blue line (down state)) converge to 
upper steady state at X T = 95 and lower state at X T = 0.0001. 



We obtain these different conditions by setting the differ- 
ent values of auto-phosphorylation parameter k 3 . For 
condition I we set k 3 = SOO^M 1 s" 1 , whereas for condition 
II we set this parameter at 0.5 fiM' 1 s" 1 . The values of all 



the other parameters are identical in both conditions 
(Table 1). For all these cases a 10 second (Ca +2 ) 4 -CaM 
pulse is used to provide the necessary stimulus to drive 
system from the basal to the up-state. 



Table 1 The numerical values of different parameters of bistable molecular loop. 



Parameter 


Description 


Value/Range 


Ai 


CaMKII degradation rate (In-active state) 


[0.0001 (sec- 1 )] 0.001 -0.00002(sec- 1 ) 


A 2 


CaMKII degradation rate (Active, phosphorylated) 


[0.0001 (sec- 1 )] 0.001 -0.00002(sec- 1 ) 


ki 


Association rate constant (Ca +2 ) 4 -CaM/CaMKII binding 


0.0011 (uM- 1 sec- 1 ) 


k2 


Dissociation rate constant CaCaM.CaMKII 


14 (sec- 1 ) 


ki 


Association rate constant CaMKII A P /CaMKII A P binding 


0.085; (uM~ 1 sec- ) 


k 2 


Dissociation rate constant CaMKII A P .CaMKII A P 


0.143 (sec" 1 ) 


k 3 


Rate constant for Phosphorylation 


500 (sec- 1 ) 


k4 


Rate constant for CaMKII A P dephosphorylation 


0.0012 (sec- 1 ) 


k 5 


Association rate constant CaMKII A P /CPEB1 binding 


0.0072;(uM- 1 sec- 1 ) 


k55 


Dissociation rate constant CaMKII A P .CPEB1 


20(sec- 1 ) 


k 6 


Rate constant for CPEB1 Phosphorylation 


0.962(sec- 1 ) 


klO 


Association rate constant (Ca +2 ) 4 -CaM/CaMKII binding 


0.001 (mM 1 sec- 1 ) 


klOlO 


Dissociation rate constant CaCaM.CaMKII 


0.8 (sec- 1 ) 


k 7 


Rate constant for CPEB1 dephosphorylation 


0.012(sec- 1 ) 


k 8 


Association rate constant PCPEB/E binding 


0.08(mM- 1 sec- 1 ) 


k88 


Dissociation rate constant PCPEB1.E 


10(sec- 1 ) 


k 9 


Synthesis Rate of New CaMKII molecule through polyadenylation 


0.08(sec- 1 ) 


Xbasal 


Basal concentration of CaMKII 


0.0001 uM 
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For condition I the dynamics from all three alternative 
models converge to same upper and lower steady states 
(X T = 95 for upper steady state and X T = 0.0001 for 
lower steady state Figure 3a. The full model is repre- 
sented by solid red line (up state) and solid blue line 
(down state). The various dashed lines present the lower 
dimensional approximations. We see that even though in 
condition I the steady state behavior is well approximated 
by the reduced systems, the dynamics are not. By setting 
derivatives to zero and replacing some dynamical vari- 
ables by functions we have produced reduced dynamics 
that are faster than the full model dynamics. These types 
of approximations can produce reasonable dynamics if 
there is a clear separation of time scales, in which case 
setting the derivatives of the fast dynamics to zero causes 
only minor differences in the dynamics. Here, although 
degradation is the slowest dynamical variable, other pro- 
cesses such as dephosphorylation can also be slow, and 
therefore there is no clear time scale separation. For con- 
dition II (Figure 3b) even the steady state behavior of the 
systems is not well approximated by the reduced 
equations. 

Next, we analyze the dependence of the steady state 
solution of this loop on different system parameters. This 
is similar to what we do in Figure 2, but here we show the 
full bifurcation diagrams resulting from a change in para- 
meters. We also compare cases where the approximations 
used are appropriate or not. First, we select the degrada- 
tion rate (A) of aCaMKII as a bifurcation parameter. Here, 
degradation represents the removal of aCaMKII either by 
general cellular degradation pathway or by diffusion from 
a certain specific location e.g., active synapses. We ana- 
lyzed the effect of degradation parameter on the character- 
istics of polyadenylation loop under the two different set 
of parameters (conditions I and II). We show results of 
this analysis in Figure 4. On the left we show the G' and F 
functions (Figure 4a,c), and on the right the complete 
bifurcation diagrams in terms of A (Figure 4b,d). Figures 4 
a, b are for condition I and 4 c, d is for condition II. 

We show the F functions for two different values of A, 
where "1" denotes X = 0.0001s" 1 and "2" denotes X = 
0.0003s" 1 . By this method we can graphically locate the 
values of X T at the fixed points. For condition I the upper 
stable steady state at X T = 95, while the lower and unstable 
steady state are at X T = 0.0001 and X T = 9.4 when degrada- 
tion rate is set at 0.0001 s" 1 (curve "1", Figure 4a). As the 
degradation rate is increased to 0.0003 s" 1 (curve "2", Fig- 
ure 4a) the new solution is located at X T = 31 (upper stable 
steady state), X T = 0 (lower stable steady state) and X T = 
9.2 (unstable steady state). Numerical bifurcation diagrams 
(Figure 4b) are developed through tracking the steady state 
behavior of all three levels of simplification with respect to 
the degradation rate of aCaMKII. The solutions contained 
in these bifurcation diagrams are tracked through a 



numerical bifurcation package [30], for the 6 and three 
dimensional models, by simply finding the cross-over in the 
one dimensional model, and by finding zeros of a polyno- 
mial in the polynomial approximation. The four bifurcation 
diagrams are nearly identical for the entire range of degra- 
dation parameter. In contrast, to these results when simula- 
tions are carried out with second set of parameters 
(condition II) the steady state solution characteristics of full 
scale model (Figure 4d dotted blue line represent the 
bifurcation diagram based on full model) does not 
match either with the three state reduced model or single 
equation model or the polynomial model (Figure 4d solid 
blue lines represent the bifurcation diagrams based on 
reduced three differential equation and a single equa- 
tion model). What these results also indicate is that the 
most significant approximation made here is in equation 1, 
and it fails when the conditions of this approximation are 
not met. 

We also analyzed the steady state solution characteristics 
of this molecular switch with respect to CPEB1 activation 
parameter k 5 . First, the simulations are implemented with 
parameters from condition I (Figure 5, a, b) and then with 
condition II (Figure 5, c, d). Here, again the steady state 
solution characteristics of the full scale model do not 
match (Figure 5d dotted blue line represent the bifurca- 
tion diagram based on full model, whereas solid blue 
lines represent the bifurcation diagrams based on 
reduced three differential equation and a single equa- 
tion model) with reduced versions when condition two is 
implemented, however, with parameters representing the 
condition one there is a complete match (Figure 5b dotted 
blue line represent the bifurcation diagram based on 
full model, whereas solid blue lines represent the bifur- 
cation diagrams based on reduced three differential 
equation and a single equation model). 

Discussion 

We have postulated that a translation-activation loop can 
form a bistable switch that could explain the mechanism 
for maintaining long term memories. Here, we explore a 
specific case of this hypothesis composed of a positive feed- 
back loop between aCaMKII and its translation regulator 
CPEB1. The possibility of a positive feedback loop is con- 
firmed by two previous experimental observations (a) phos- 
phorylation of CPEB1 regulates the synthesis of aCaMKII 
molecules through polyadenylation of aCaMKII-mRNA 
[20] (b) aCaMKII phosphorylates the CPEB1 molecule 
[25,26]. Thus, both the aCaMKII and CPEB1 interact 
through a closed loop. Based on the elementary enzyme 
kinetics (Method A) the dynamics of this loop can be char- 
acterized through a high dimensional system of ordinary 
differential equations. One can study the characteristics of 
such a system by numerically integrating these differential 
equations and with different initial conditions and kinetic 
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Figure 4 Steady state solution characteristics of molecular loop with respect to parameter A. This contain an analytical solution (a, c) and 
numerical/analytical bifurcation diagrams (b, d). Approximate analytical solution is developed through graphing equation 6, whereas the 
numerical/analytical bifurcation diagrams are developed through tracking the steady state behavior of all three models with respect to the 
degradation rate. Two different set of parameters are compared. For first set of parameters the approximate solution and numerical/analytical 
bifurcation diagrams from three alternative models converge to same upper, stable steady state solution branch (a, b). For the second set of 
parameters the approximate solution and numerical/analytical bifurcation diagrams from three alternative models does not converge to same 
upper, stable steady state and solution branch(c, d). The approximate analytical solution for first set of parameters (a) is located at two different 
degradation rates (The solid blue line#"1" represent A = 0.0001s 1 and other solid blue line#"2" represent A = 0.0003 s 1 ). This method 
locates the upper stable steady state at X T = 95, while the lower and unstable steady state are at X T = 0 and X T = 9.4 when degradation rate is 
set at 0.0001 s" 1 (curve 1). As degradation rate is increased to 0.0003 s" 1 (curve 2) the new solution is located at X T = 31 (upper stable steady 
state), X T = 0.0001 (lower stable steady state) and X T = 9.2 (unstable steady state). The numerical and analytical bifurcation diagrams for first set 
of parameters (b) are developed through tracking the steady state behavior of all three models. Three bifurcation diagrams exactly match with 
each other for the entire range of degradation parameter (dotted blue line represent the bifurcation diagram based on full model, 
whereas solid blue lines represent the bifurcation diagrams based on reduced three differential equation and a single equation 
model). The approximate analytical solution for second set of parameters (c) is also located at two different degradation rates (The solid blue 
curve "1" represent A = 0.0001s" 1 , and other solid blue curve "2" represent A = 0.0003 s 1 ). This method locates the upper stable steady 
state at X T = 95, while the lower and unstable steady state are at X T = 0.0001 and X T = 9.4 when degradation rate is set at 0.0001s- 1 (curve 1). 
As degradation rate is increased to 0.0003 s 1 (curve 2) the new solution is located at X T = 31 (upper stable steady state), X T = 0.0001 (lower 
stable steady state) and X T = 9.2 (unstable steady state). The numerical and analytical bifurcation diagrams for second set of parameters (d) are 
developed. The bifurcation diagram from full model does not match with the bifurcation diagrams from reduced and single equation model. 

>x . J 



rate parameters, and extensively search for the possibility 
of a two distinct stable steady states. This procedure itself 
is very tedious, time consuming and inaccurate. 

Here, we develop a systematic approach to reduce the 
dimensionality of a high dimensional ODE based model 
of aCaMKII-CPEBl molecular loop. Our simplification 
strategy is devised such that the essential features of full 
scale model are preserved. This simplification process is 
based on introducing the two key approximations. 1. The 
synthesis-degradation curve (equation 1) and 2. The fit- 
ting function curve (equation 7). Through these approxi- 
mations the full scale model is reduced to three lower 
dimensional models, (a) A three dimensional model com- 
prising of equations 6, B-6, and B12. (b) A one dimen- 
sional model (equation 6), (c) A one dimensional fourth 
order polynomial based model (B-18). We can extract the 
fixed point solution characteristics of this system from all 



three reduced models however, only reduced model (a) 
and (b) can provide the dynamical characteristics of this 
system. Since, during the simplification process it is pos- 
sible to loose some features therefore, it would be logical 
to compare the performance of these three reduced mod- 
els with a full scale model of this loop. 

Our results show that the key to performance of three 
reduced models in locating fixed points of the system is 
the approximation introduced through equation 1. Since, 
in this approximation the CaMKII bound in biochemical 
complexes also degrade along with free CaMKII, which is 
in contrast to full scale model. Thus, for parametric condi- 
tions where the amount of bound CaMKII is negligible the 
three reduced models yield matching results to full scale 
model. However, for conditions where there is a non-neg- 
ligible amount of bound CaMKII the results of three mod- 
els are not matching with full model. Interestingly, 



Aslam and Shouval BMC Systems Biology 2012, 6:12 
http://www.biomedcentral.eom/1 752-0509/6/1 2 



Page 11 of 13 






100 




80 




60 


|| 


40 


Coi 
Tut 


20 





d 



^ | 300 

! 1 200 

§ P 

S I ioo 



.'I I 



Figure 5 Steady state solution characteristics of molecular loop with respect to CPEB1 activation parameter k 5 The parametric steady state 
solution characteristics are developed through an approximate analytical solution (a, c) and numerical/analytical bifurcation diagrams (b, d). The 
approximate analytical solution is developed through graphing equation 6, whereas the numerical/analytical bifurcation diagrams are developed 
through tracking the steady state behavior of all three models with respect to k 5 . Two different set of parameters are compared. For first set of 
parameters the approximate solution and numerical/analytical bifurcation diagrams from three alternative models converge to same upper, stable 
steady state solution branch (a, b). For the second set of parameters the approximate solution and numerical/analytical bifurcation diagrams from 
three alternative models does not converge to same upper, stable steady state and solution branch (c, d). The approximate analytical solution for first 
set of parameters (a) is located at two different values of rate constant k5 (The dotted red curve "1" represent k5 = 0.0072 jj\\A'\s'\ second dotted 
red curve "2" represent k5 = 0.0001 ^/M~ 1 .s" 1 ). This method locates the upper stable steady state at X T = 95, while the lower and unstable steady 
state are at X T = 0.0001 and X T = 9.4 when k 5 is set at 0.0072^/M" 1 .s" 1 (curve "1"). As k 5 is decreased to 0.0001 new solution is located at X T = 70 
(upper stable steady state), X T = 0.0001 (lower stable steady state) and X T = 16 (unstable steady state). The numerical and analytical bifurcation 
diagrams for first set of parameters (b) are developed through tracking the steady state behavior of all three models with respect to k5. The three 
bifurcation diagrams are exactly matching with each other for the entire range of activation parameter (dotted blue line represent the bifurcation 
diagram based on full model, whereas solid blue lines represent the bifurcation diagrams based on reduced three differential equation 
and a single equation model). The approximate analytical solution for second set of parameters (c) is also located at two different values of rate 
constant k 5 (The dotted red curve "1" represent k 5 = 0.0072 iJ.hA~\s \ second dotted red curve "2" represent k5 = 0.0001 (i/M'V 1 ). The numerical 
and analytical bifurcation diagrams for second set of parameters (d) are developed through tracking the steady state behavior of all three models with 
respect to k5. Full model bifurcation diagram does not match with reduced and single equation model. 



however the results of these three reduced models are 
matching with each other. Here, this phenomenon of 
CaMKII trapping in biochemical complexes is simply 
modeled through an autophosphorylation parameter k 3 . 
For the k 3 set at lower values the amount of CaMKII 
bound is large and thus leading to disagreements between 
the results from full scale model and reduced models. 

The analysis of this paper also provides some intuition 
on the complex process of polyadenylation based transla- 
tion. Here, through a systematic simplification process 
we develop a single variable dynamical equation (eq. 6) 
which shows that the change in total CaMKII concentra- 
tion is due to a balance between protein synthesis and 
degradation. At steady state this equation is composed of 
two terms i.e., synthesis and degradation term. Thus, the 
polyadenylation based switching is essentially due to a 
balance between the synthesis of new proteins and degra- 
dation of old ones. If the amount of synthesis is at much 



faster rate then the degradation rate the balance will shift 
towards single state high protein concentration, similarly 
if the degradation rate is much faster then new protein 
synthesis the balance will tip towards single lower state 
of switch. So for reversible switching a balance should be 
maintained between protein synthesis and degradation. 
This equation also provides an easy method to analyze 
the effect of other system variables on switching charac- 
teristics. For example one parameter which could criti- 
cally affect the polyadenylation switching is the CPEB1 
activation rate through CaMKII i.e., k 5 . For larger value 
of this parameter the new synthesis rate will increase and 
at lower value the new synthesis rate will decrease. 

In this paper, we simplify the high dimensional ODE 
model of aCaMKII-CPEBl loop in such a way that these 
simplifications preserve the characteristics properties of 
the full scale ODE model. By reducing the high dimen- 
sional model to the minimal plausible scenario, we were 
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able to obtain a single algebraic expression that charac- 
terizes the fixed points of these dynamics. We then used 
these algebraic equations to develop the analytical bifur- 
cation diagrams by perturbing the various parameters (X 
and k5). Using this analytical expression we explain how 
the loss of aCaMKII in synapses due to protein degrada- 
tion is balanced by synthesis of new aCaMKII molecules. 
This analysis explains how the pair of two molecules in 
active and inactive forms and a synthesis based positive 
feedback loop can lead to a bistable switch. On the basis 
of these bifurcation diagrams we show that even the sim- 
plified version of a polyadenylation based model of 
aCaMKII can exhibits an up and down state and by per- 
turbing these parameters the system can toggle between 
an up and a down state. The activity induced increase in 
total amount of aCaMKII can be maintained for long 
period of time due to a balance between the new synth- 
esis and degradation of aCaMKII molecule. As we set 
equation 1 equals to zero the resulting equation 2 
describe the steady state value of total amount of aCaM- 
KII. It also shows that total concentration of CaMKII is 
directly proportional to the amount of CPEB1 phos- 
phorylated and is inversely proportional to the degrada- 
tion rate. Therefore, when the fraction of CPEB1 
phosphorylated is high the amount of CaMKII generated 
through polyadenylation increases, which will out balance 
the amount of aCaMKII removed through degradation 
and there is a net up-regulation of aCaMKII concentra- 
tion. Similarly this equation also explains that if degrada- 
tion rate is too fast the compensation provided by new 
synthesis of aCaMKII molecules will not be enough to 
balance the loss of CaMKII due to degradation. Our sim- 
plification captures the essential properties of protein 
translation through polyadenylation. It shows that how 
the activity induced signal converts the inactive form of 
aCaMKII into active, which is further amplified through 
an auto-phosphorylation loop. The active and phos- 
phorylated aCaMKII in turn drives the synthesis of a 
new inactive aCaMKII molecule, through phosphorylat- 
ing the CPEB1. Here, all these steps are analytically pro- 
ven as shown by the equations (B1-B9). 

Apart, from the approximation introduced through 
equation 1 another approximation is introduced into this 
analysis when we tried to develop a single algebraic equa- 
tion 6 based version of this system. In equation 6 we 
developed an approximate function XP = g ( XT). In 
order, to develop a completely accurate solution this func- 
tion should be analytically developed, however, for this 
system it was not possible to invert the highly non-linear 
system of XT as a function of XP. Here, we developed the 
approximate function XP = g (XT) through a fitting func- 
tion routine in matlab. The approximate fitting function 
proved to be very sufficient good in approximating both 
the upper, lower stable steady state solution and unstable 



steady state solution., however, it could not accurately 
approximate the lower stable steady state solution. 

The numerical values of different parameters of this 
molecular loop are described in table 1. Many of these 
parameters are extracted from previous observations based 
on experimental and simulation work [16,19,20,25-34]. For 
example parameters like aCaMKII degradation rates, cal- 
cium/calmodulin binding and un-binding rate, rate para- 
meters for aCaMKII auto-phosphorylation loop and rate 
parameters describing the new synthesis of aCaMKII pro- 
tein are taken from previous experimental and simulation 
based observations [27,28,31-33], and [34]. Some other rate 
parameters such as activation of CPEB1 through phos- 
phorylated aCaMKII are obtained through parameters 
scaling and matching the observed experimental dynamics. 
The parameters of table 1 and bifurcation diagrams of 
Figure 4 and 5 indicate that this molecular system is very 
robust. For example the k 5 parameter (describing aCaMKII 
mediated activation of CPEB1) exhibits the bistable charac- 
ter over a range of three orders of magnitude. This means 
that even a small activation of CPEB1 through aCaMKII 
will induce the bistable character to this loop and large 
strength of aCaMKII activity inhibition is required to 
reverse the up-state to down. This leads to interesting 
observation where many of aCaMKII inhibitors are not 
able to reverse the established L-LTP when applied during 
the maintenance phase of memory formation [33]. The 
reduced model of this work and full-scale model [30] can 
make these predictions. 

Conclusions 

This study provides a systematic method to simplify, 
approximate and analyze a molecular model of polyadeny- 
lation loop. This model of polyadenylation loop is based 
upon de-novo synthesis of new proteins & their activation. 
The polyadenylation loop operates through a positive feed- 
back loop between a protein and its translation factor. 
This work shows how to extract low dimensional systems 
that can be used to obtain analytical solutions for the fixed 
points of the system and to describe the dynamics of the 
system. The methods used here have general applicability 
to the formulation and analysis of many molecular 
networks. 
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